home *** CD-ROM | disk | FTP | other *** search
/ MACD 5 / MACD 5.bin / workbench / libs / unixlib.lha / unix / src / hypot.c < prev    next >
C/C++ Source or Header  |  1995-12-30  |  7KB  |  252 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.  
  34. /* 
  35.  * IEEE 754 recommended functions:
  36.  * (a) copysign(x,y) 
  37.  *              returns x with the sign of y. 
  38.  * (b) scalb(x,N) 
  39.  *              returns  x * (2**N), for integer values N.
  40.  * (c) logb(x) 
  41.  *              returns the unbiased exponent of x, a signed integer in 
  42.  *              double precision, except that logb(0) is -INF, logb(INF) 
  43.  *              is +INF, and logb(NAN) is that NAN.
  44.  * (d) finite(x) 
  45.  *              returns the value TRUE if -INF < x < +INF and returns 
  46.  *              FALSE otherwise.
  47.  *
  48.  *
  49.  * CODED IN C BY K.C. NG, 11/25/84;
  50.  * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
  51.  */
  52.  
  53. #include <math.h>
  54.  
  55. static const unsigned short msign = 0x7fff, mexp = 0x7ff0;
  56. static const short prep1 = 54, gap = 4, bias = 1023;
  57. static const double novf = 1.7E308, nunf = 3.0E-308, zero = 0.0;
  58.  
  59. double
  60. scalb(double x, int N)
  61. {
  62.         int k;
  63.  
  64.         unsigned short *px = (unsigned short *) &x;
  65.  
  66.         if (x == zero)
  67.         return(x); 
  68.  
  69.         if ((k = *px & mexp) != mexp) {
  70.             if (N < -2100)
  71.             return(nunf*nunf);
  72.         else if (N > 2100)
  73.             return(novf+novf);
  74.             if (k == 0) {
  75.                  x *= scalb(1.0, (int)prep1);
  76.          N -= prep1;
  77.          return(scalb(x, N));
  78.         }
  79.  
  80.             if((k = (k >> gap) + N) > 0) {
  81.                 if (k < (mexp>>gap)) {
  82.             *px = (*px & ~mexp) | (k << gap);
  83.                 } else {
  84.             x = novf+novf;               /* overflow */
  85.         }
  86.             } else {
  87.                 if (k > -prep1) {                /* gradual underflow */
  88.                     *px = (*px & ~mexp) | (short)(1 << gap);
  89.             x *= scalb(1.0, k-1);
  90.         } else {
  91.                     return(nunf*nunf);
  92.         }
  93.         }
  94.         }
  95.         return(x);
  96. }
  97.  
  98.  
  99. double
  100. copysign(double x, double y)
  101. {
  102.         unsigned short  *px = (unsigned short *) &x,
  103.                         *py = (unsigned short *) &y;
  104.  
  105.         *px = (*px & msign) | (*py & ~msign);
  106.         return(x);
  107. }
  108.  
  109. double
  110. logb(double x)
  111. {
  112.     short *px = (short *)&x, k;
  113.  
  114.     if ((k = *px & mexp) != mexp) {
  115.         if (k != 0) {
  116.             return((double)((k >> gap) - bias));
  117.         } else if (x != zero) {
  118.             return(-1022.0);
  119.         } else {
  120.             return(-(1.0/zero));
  121.         }
  122.     } else if (x != x) {
  123.         return(x);
  124.     } else {
  125.         *px &= msign;
  126.         return(x);
  127.     }
  128. }
  129.  
  130. int
  131. finite(double x)
  132. {
  133.         return( (*((short *) &x ) & mexp ) != mexp );
  134. }
  135.  
  136. /* HYPOT(X,Y)
  137.  * RETURN THE SQUARE ROOT OF X^2 + Y^2  WHERE Z=X+iY
  138.  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  139.  * CODED IN C BY K.C. NG, 11/28/84;
  140.  * REVISED BY K.C. NG, 7/12/85.
  141.  *
  142.  * Required system supported functions :
  143.  *    copysign(x,y)
  144.  *    finite(x)
  145.  *    scalb(x,N)
  146.  *    sqrt(x)
  147.  *
  148.  * Method :
  149.  *    1. replace x by |x| and y by |y|, and swap x and
  150.  *       y if y > x (hence x is never smaller than y).
  151.  *    2. Hypot(x,y) is computed by:
  152.  *       Case I, x/y > 2
  153.  *        
  154.  *                       y
  155.  *        hypot = x + -----------------------------
  156.  *                         2
  157.  *                sqrt ( 1 + [x/y]  )  +  x/y
  158.  *
  159.  *       Case II, x/y <= 2 
  160.  *                                   y
  161.  *        hypot = x + --------------------------------------------------
  162.  *                                       2 
  163.  *                                 [x/y]   -  2
  164.  *               (sqrt(2)+1) + (x-y)/y + -----------------------------
  165.  *                                       2
  166.  *                              sqrt ( 1 + [x/y]  )  + sqrt(2)
  167.  *
  168.  *
  169.  *
  170.  * Special cases:
  171.  *    hypot(x,y) is INF if x or y is +INF or -INF; else
  172.  *    hypot(x,y) is NAN if x or y is NAN.
  173.  *
  174.  *    hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units
  175.  *    in the last place). See Kahan's "Interval Arithmetic Options in the
  176.  *    Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
  177.  *    1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
  178.  *    code follows in comments.) In a test run with 500,000 random arguments
  179.  *    on a VAX, the maximum observed error was .959 ulps.
  180.  *
  181.  * Constants:
  182.  * The hexadecimal values are the intended ones for the following constants.
  183.  * The decimal values may be used, provided that the compiler will convert
  184.  * from decimal to binary accurately enough to produce the hexadecimal values
  185.  * shown.
  186.  */
  187.  
  188. const static double r2p1hi = 2.4142135623730949234E0;
  189. const static double r2p1lo = 1.2537167179050217666E-16;
  190. const static double sqrt2 = 1.4142135623730951455E0;
  191.  
  192. double
  193. hypot(double x, double y)
  194. {
  195.     static const double zero = 0.0, one = 1.0,
  196.             small = 1.0E-18;    /* fl(1+small)==1 */
  197.     static const ibig = 30;    /* fl(1+2**(2*ibig))==1 */
  198.     double t, r;
  199.     int exp;
  200.  
  201.     if (finite(x)) {
  202.         if (finite(y)) {    
  203.         x = copysign(x, one);
  204.         y = copysign(y, one);
  205.         if (y > x) {
  206.             t=x;
  207.             x=y;
  208.             y=t;
  209.         }
  210.         if (x == zero)
  211.             return(zero);
  212.         if (y == zero)
  213.             return(x);
  214.         exp = logb(x);
  215.         if (exp-(int)logb(y) > ibig ) {
  216.             /* raise inexact flag and return |x| */
  217.             one+small;
  218.             return(x);
  219.         }
  220.  
  221.             /* start computing sqrt(x^2 + y^2) */
  222.         r = x-y;
  223.         if (r > y) {        /* x/y > 2 */
  224.             r = x/y;
  225.             r = r+sqrt(one+r*r);
  226.         } else {        /* 1 <= x/y <= 2 */
  227.             r /= y;
  228.             t = r*(r+2.0);
  229.             r += t/(sqrt2+sqrt(2.0+t));
  230.             r += r2p1lo;
  231.             r += r2p1hi;
  232.         }
  233.  
  234.         r = y/r;
  235.         return(x+r);
  236.  
  237.         } else if (y == y) {    /* y is +-INF */
  238.         return(copysign(y, one));
  239.         } else {
  240.         return(y);        /* y is NaN and x is finite */
  241.         }
  242.     } else if (x == x) {        /* x is +-INF */
  243.         return (copysign(x, one));
  244.     } else if (finite(y)) {
  245.         return(x);            /* x is NaN, y is finite */
  246.     } else if (y != y) {
  247.         return(y);            /* x and y is NaN */
  248.     } else {
  249.         return(copysign(y, one));    /* y is INF */
  250.     }
  251. }
  252.